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CONTINUOUS ANALOG OF THE GAUSS-NEWTON METHOD 
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A Continuous Analog of discrete Gauss-Newton Method (CAGNM) for numerical solu- 
tion of nonlinear problems is suggested. In order to avoid the ill-posed inversion of the 
Frechet derivative operator some regularization function is introduced. For the CAGNM 
a convergence theorem is proved. The proposed method is illustrated by a numerical 
example in which a nonlinear inverse problem of gravimetry is considered. Based on 
the results of the numerical experiments practical recommendations for the choice of the 
regularization function are given. 
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1. Introduction 



Let Hi and H2 be real Hilbert spaces and (p : H\ — > H2 a nonlinear operator. 
Let us consider the equation: 

<p(x)=0. (1.1) 
We assume that the following condition on <p holds. 



Condition A: Problem (LI) has a solution x, not necessarily unique. 

In the well-known Newton's method (!) one constructs a sequence {x n } for 
n = 0, 1, . . . which converges to a solution (in general non unique) of equation 



(1.1). The first term Xq is an initial approximation point and the other terms are 



constructed by means of the following iterative process: 

x k+1 =x k -(p'{x k y L (p{x k ), (1.2) 

where ip'(x) is the Frechet derivative of the operator ip. Recall that <p'(x) is a 
linear operator from H\ to 7?2- The usual necessary condition for the realization 
of the Newton method is the bounded invertibility of (p'(x k ), that is, the existence 



*E-mail: airapet@math.ksu.edu 
^E-mail: ramm@math.ksu.edu 
t E-mail: smirn@math.ksu.edu 



of a bounded linear operator [(//(a^)] -1 for all k. Actually in order to provide the 
convergence of Newton iterations one needs bounded invertibility of ip' in a ball 
B(x,R) :— {x : x € Hi, \\x — x\\ < R}. However this condition does not hold in 
many important applications. In order to avoid this restriction several modifications 
of the Newton method have been developed. In this paper we consider the Gauss- 
Newton procedure for equation ( |1.1| ) (see e.g. (□)): 

x k+ i = x k - [ip'*(x k )Lp'(x k )]~ 1 ip'*{x k )ip(xk), Xo = X . (1.3) 

If the operator ip'*(x)(p'(x) is not boundedly invertible one needs some regularization 
procedure. In order to construct such a procedure one can introduce a sequence of 



positive numbers a k , a k — > 0, and replace iterative method (L3) by the following 
li): 



one 



x k+ i =x k - [ip'*(x k )tp'(x k ) + a k I] 1 [tp'*(x k )(p(x k ) + a k {x k - x )], (1.4) 

where / is the identity operator. 

The methods constructed above can be also considered as discrete analogs of 
some continuous methods (called sometimes continuation methods). In (til) the 
following Cauchy problem has been considered as a continuous analog of (fT^): 

dx 

x(t) = - l p'(x(t)Mx(t)), x(0) = x , ±(t):= — . (1.5) 



A solution to problem (1.1) can be obtained as limit of the function x(t) for t — > oo. 



If one solves this Cauchy problem by means of Euler's method with a stepsizc r = 1 



one gets (1.2) with x k = x{k). Continuous analogs of iterative methods have several 
advantages over the discrete ones. Convergence theorems for continuous methods 
usually can be obtained easier. If a convergence theorem is proved for a continuous 
method, that is, for the Cauchy problem for a differential equation, for instance 



(1.5), one can construct various finite difference schemes for the solution of this 



Cauchy problem. These difference schemes give discrete methods for the solution 
of equation (1.1). For instance the methods of Euler and Runge-Kutta can be used. 
More detailed information about the. applications and modifications of continuous 
Newton methods can be found in (uOEJ). 

The aim of this paper is to construct a continuous analog of iterative scheme 



(1.4), to prove a convergence theorem for this continuous analog of (1.4), and to 
test the method numerically by applying it to a practically interesting nonlinear 
inverse problem of gravimetry. 

The paper is organized as follows. In section 2 a continuous analog of method 



(1.4) is described and a convergence theorem for this method is formulated. In 
section 3 this convergence theorem is proved. In section 4 an inverse gravimetry 
problem is considered and the proposed method is numerically tested. In our nu- 
merical experiments comparison of different regularization functions is done. Based 
on the results of the numerical experiments some recommendations are given for 
the choice of the regularization function. 



2. Continuous Gauss-Newton Method and Convergence Theorem 



In order to describe convergence rates we introduce the following 

Definition 2.1. A positive function a(t) G C 1 [0,oo) is said to be a convergence 
rate function if a(t) decreases monotonically to zero as t — > oo, a(t) G C 1 [0, oo) and 
lna(i) is concave, that is, a(t)/a(t) is monotonically increasing. 

Remark 2.2. The number a(0) can be chosen sufficiently large and simultaneously 
the number |a(0) / a(0)| can be sufficiently small. Here and below the over dot 
denotes the derivative with respect to time x :— dx/dt. For example, one can 
choose a(t) — b/(t + a), where a and b are positive constants such that a and b/a 
are sufficiently large. 

A continuous analog of iterative process (|l.4|) is the following Cauchy problem: 
x(t) = -[<(/• (x(t))<f/{x(t)) + a{i)I]- x W*{x{t))y{x{t)) + a(t)(x(t) - x )], (2.1) 

x(0) = x a . 

Denote by Ran(L) the range of the linear operator L. The convergence of the con- 
tinuous analog of Gauss-Newton method (CAGNM) is established by the following 
theorem, in which (and throughout this paper) the norms ||y'(a;)|| and ||y"(x)|| are 
the norms of linear and bilinear operators from Hi to H2 and from Hi x Hi to H2 
respectively. 

Theorem 2.3. Let a(t) be a convergence rate function. Assume that there exists 
a positive number R for which Condition A and the following conditions hold: 

(i) The Frechet derivatives f'(x) and f"(x) exist in the ball B(x,R) and satisfy 
the following inequalities: 



where 



\\v'(x)\\ <N 1} \\<p"(x}\\ < N2 \JxeB(x,R), (2.2) 



a(0) (l-2NiN 2 \\v\\ + ^-)<R. 



NiN 2 V a(0) 



(ii) 



x G B(x , R)n[x + Ran((p'*{x)(p'(x ))}. (2.3) 
(iii) For some v, such that x — xq = ip'* (x)ip' (x)v, the following inequalities hold: 

l-2NiN 2 \\v\\ + ^->0, (2.4) 

l_2JV 1 jyr 2 |H| + ^ > ) - 2NiN 2 \\v\\ > 0. (2.5) 
a{0)J 



Then the following conclusions hold: 



(i) The solution x — x(t) of problem (2.1) exists, and x(t) £ B(x, R) for t £ [0, oo), 

(ii) \\x(t) - x\\ = 0{a(t)) for t -> oo. 

Remark 2.4. Condition (ii) in Theorem 2.3 gives some restriction on the choice of 
an initial approximation point. It is not easy to verify this condition algorithmically. 



However some kind of this condition is necessary if one works with equation (1.1) 
with the operator ip'*(x)tp' (x), which is not boundedly invertible. If the operator 
ip'*<p' is injective but is not boundedly invertible, then the image of the linear 
selfadjoint operator (f'*(f' is dense in B(x, R) and consequently the set of the suitable 
initial approximation points satisfying condition (ii) is also dense in B(x, R). As our 
numerical results show (see section 4) the proposed method is practically efficient. 

3. Proof of Theorem 2.3 



The main part of the proof is to show that the solution to problem (2.1) does 
not leave the ball B(x,R) (Lemma 3.3). In order to prove it, let us assume that 
there exists such a point t% £ [0, oo) that x(t) intersects the boundary of B(x,R) 
for the first time at t = t%. Hence x(t) belongs to the interior of the B(x,R) for 
t G [0,ti) and ||x(fi) — x\ \ = R. Let us introduce an auxiliary function 

w (t) := \\ x (t) -x\\/a{t). (3.1) 

First, in Lemma 3.1, we derive a nonlinear differential inequality for w(t). From 
this differential inequality we get the estimate which shows that for all t £ [0, t\] 



the points of the integral curve of problem (2.1) belong to the interior of the ball 
B(x,R). This contradiction proves that the integral curve of the solution does not 
leave the above ball, and consequently problem ( |2.l| ) has the global solution for 
t £ [0, oo). Also we show the boundedness of the w(t) and this implies, by formula 



(3.1), strong convergence of x(t) to x for t — * oo. 



Lemma 3.1. If the assumptions of Theorem 2.3 hold then the differential inequal- 
ity 

^ <C lW 2 -C 2 w + C 3: (3.2) 
at 

is valid for t £ [0,ii], where 

Ci = ^A Q l = l-2N 1 N2\\v\ + =^-, C 3 = v. (3.3) 



Proof. The Gateaux derivative ^1,^2) is a bilinear operator such that 

¥>'(*+£i)k-^(*)6 := and 1 1 ^7 1 1 - 1 1 1 1 " 1 -> for ^ - 0, \M > 0. 



Let us define operators K, G : B(x, R) X H \ X H \ — * H 2 by the formulas: 



and 



1 1 

K(x,£i,&) = J J(f/'(x + 8tZi,Zi,&)tdtd8 


1 



(3.4) 



(3.5) 



Then from (2.2) we get 



No 



\K{x,^,i 2 )\ < -f\M ■ H6II, \G(x,^ 2 )\ < N 2 \M ■ H6I 
The following formulas will be used: 

<p(x) — ip(x) = ip'{x)(x — x) + K(x, x — x, x — x), 

and 

(<//(£) - <p'(x))£ = G(x,x-x,g), 



(3.6) 



(3.7) 



(3.8) 



where K and G are defined by (3.4) and (p.q) respectively. Let us derive formulas 



(3.7) and (fej. One has 



(p(x) — <p(x) — I —ip(x + t(x — x))dt = I ip'(x + t(x — x))(x — x)dt 



<p'(x)(x — x) + \<f'{x + t(x — x)) — <p'(x)](x — x)dt 



1 1 

(p'(x)(x — x) + J tdt J ds[(p"(x + st(x — x))](x — x)(x — x) 




and 



[<p'(x) - <p'(x)]£= j —<p'(x + t(x-x))£dt = I <p"(x + t(x-x),x-x,£)dt 



Since x solves (1.1), one can rewrite equation (2.1) as 
dx 

— = -[ip'*{x)ip'(x)+al]- 1 [<p'*{x)(<p(x) -<p(x))+a(x-x )]. 



From the condition (iii) of Theorem 2.3 and from (|3.7| ) we get 

dx 

— = -[<p'*(x)tp'(x) + oJ] -1 [-<,£>'* (x)(<p' (x)(x -x) + K(x,x- x,x- x))+ 

(XL 



a(x — x) + aip'* (x)ip' (x)v], 

and therefore 

— = — (x — x) — [if'* (x)(p' (x) + oel]^ 1 [— ip'*(x)K(x, x — x, x — x) + a<f'*(x)(f'(x)v+ 

CLT 

a{if'*{x)if'{x)-ip'*{x)f'{x))v]. 
Since x does not depend on t, it follows from ( |3.8| ) that 

X) = ^j- = -(x-x)- [f'* {x)f'{x) + aI]- 1 [~f'*(x)K(x, x-x,x-x) 

+aip'* {x)if' {x)v + a((f'* (x) — if'* (x))if'(x)v + aip'* {x)G(x, x — x, v)]. 
Now let us derive an inequality for — x\\ 2 . One has 

— \\x — x\\ 2 = — 2| |x — x\ | 2 + 2([(p'* (x)ip' (x) + a I] ~~ 1 [if'* (x)K(x, x — x,x — x)],x — x) — 

2a{[<f'* ix)if' (x) + al]^ 1 if'* (x)f (x)v , x — x) + 2a(<f'(x)v, G(x, x — x, [f* {x)f' {x) 

+al]~ 1 (x - x))) + 2a([f'*(x)f'(x) + al]~ V* {x)G(x, x - x,v),x- x). 

Since the operator ip'*(x)<f' (x) is selfadjoint and nonnegative we have the following 
spectral representation: 

W*{x)f'{x) + al]~ l y'*{x)<f'{x) = / t— d**, 

Jo * + a 

where E\ is the resolution of the identity of the selfadjoint operator f'*{x)f'{x). 
Since < A/(A + a) < 1 for a > and A > 0, it follows that 

\\[f'*(x)f'(x)+a(t)ir 1 f'*(x)f'(x)\\ < 1. (3.9) 

Also one has the following estimate: 

\\[f'*(x)f'(x)+a(t)I]- 1 \\<l/a(t). (3.10) 

From (p|),(p^|) and |[|) one gets the following differential inequality for A(t) :— 
\\x(t)-x[[: 

A < -A+ ^lA 2 + a\\v\\ +2JViJV 2 |M|A 
2a 

In order to finish the proof of the Lemma 3.1, we derive from the last inequality the 
inequality for w(t) by taking into account that a(t)/a(t) is monotonically increasing 
function. □ 

The following lemma is a simple corollary of the more general results established 
in ft. 



Lemma 3.2. Let f(t, u) be a continuous function on [0, T] x (—00, +00) such that 
the Cauchy problem 

ii = f(t,u(t)), u(0) = u o (3.11) 

is uniquely solvable on [0, T] and v(t) be a differentiable function defined on [0, T] 
and satisfies the conditions 

v<f(t,v(t)), te[0,T], v(0)=v o . (3.12) 

If vq < uo then 

v{t) < u(t) for t € [0,T]. 



It follows from inequality (21) that c = \J C\ — 4C1C3 > for constants Ci, 
C 2 , C3 defined in (3.3). Let u\ and u 2 be correspondingly the smaller and the 
larger roots of the equation C\u 2 — C^u + C3 = 0. For uq satisfying the inequality 
U2 < uq < C2/(2Ci) the solution of the Cauchy problem 



C 1 u 2 -C 2 u + C 3 , u(0)=u . 



(3.13) 



is given by the formula: 



u-u 2 



u — Ul 



u 2 - u gCt 



u - Ul 

Let us show that u(t) is defined for all t £ [0,oo). Indeed, uq £ (^1,^2), so for 
sufficiently small t one has 



u 2 - Ui 

U2~uq r .t 1 



(3.14) 



Thus ui < u(t) < u(0) < u 2 . This means that u(t) does not leave the interval 
(ui,u 2 ) for all t £ [0, 00) and u(t) is well defined for all t £ [0, 00). 

From condition (i) of Theorem 2.3 one obtains the following estimate: 



w(0) 



x -x\\ < l_-2N 1 N 2 \\v\\ + l$ 



a(0) iViiV 2 
Therefore from Lemmas 3.1 and 3.2 it follows that 



2Ci 



\\x{ti) — x\\ u 2 ~ ui C 2 

a( h ) ^ Wl+ ^ e - 1+ l <U0< 2Cr- 

v / 1/.n — 11.1 



(3.15) 



Thus 



This contradicts the assumption | |a;(ii) — x \ \ = R. So the following lemma is proved. 



Lemma 3.3. If the assumptions of Theorem 2.3 hold and for an arbitrary positive 
T the solution of the problem (^J) exists on the interval [0, T], then the integral 



curve of the solution of (2T) lies in the interior of the ball B(x,R) for all t from 
the interval [0,T]. 



Now let us show that there exists the unique solution of (|2.l|) on [0, oo)provided 



that xq satisfies conditions ii) and iii) of Theorem 2.3. The Cauchy problem (2.1) 
is equivalent to the integral equation 

t 

x(t)=x + J F(s,x(s))ds, (3.16) 
o 

where 

F(s,x(s)) := -[ip'*(x(s))<p'(x(s)) + q(s)/] -1 [^'* (x(s))<p(x(s)) + a(s)(x(s) - x )}. 

Let us fix an arbitrary large positive number T and use the successive approximation 
method to solve equation ( 3.1 6| ) on [0, T]: 



t 

x n+1 {t) = x + J F(s,x n (s))ds, x (t)=x , te[Q,T]. (3.17) 
o 

Since tp'(x) and tp"(x) are assumed to be bounded in B(x,R), see (|2.2D, and a(t) 
is positive on [0,T], for every t e [0, T] the function F(t,x) has bounded Frechet 
derivative with respect to x in B(x,R). So one has: 

\\F{t 7 x 1 )-F(t,x 2 )\\ < K{T)\\x x -x 2 \\ 
for all t £ [0, T] and xi, x 2 belong to B(x, R). Thus, one easily gets the estimate 

rpn 

\\x n+1 (s) - x n (s)\\ < \\F(x )\\K n (T) — 

valid on the maximal subinterval [0,Xi] = {t : t G [0, T] and x(t) £ B(x,R)}. 
Therefore iterative process fl3.17| ) converges uniformly and determines the unique 



solution of equation ( 3.1£) on [0, T{\. If Ti < T, it follows from the maximality of the 



subinterval [0, Ti] that x{T{) is a boundary point of B{x, R). But this contradicts to 



Lemma 3.3. So the solution of the problem (2T) exists and belongs to the interior 
of the ball B(x, R) on every interval [0, T] and consequently on [0, oo). 

To finish the proof of Theorem 2.3 it is sufficient to note that equation ( 3.1 5| ) 
implies estimate 



\x(t)-x\\<^-a(t) 



for all t e [0,oc). □ 



4. Numerical Results 



To test numerically the method described above, we chose the inverse gravime- 
try problem (§). The goal of the numerical test is to illustrate the choice of the 
regularization function a(t) and to compare two methods of solving the Cauchy 
problem (|]l]): the Euler method, which corresponds to the iterative scheme ( |l.4| ), 
and the Runge-Kutta method. 

Let the sources of a gravitational field with a constant density p be distributed 
in the domain 

D = {-l<t<l, -H < z < -H + x(t)}, 

where x(t) is an interface between two media, I and H are parameters of the domain. 
The potential V of such a field is given by the double integral: 

V(t,z) = ^- [ [pin — 1 dS = 



I -H+x(s) 

P_ 

47T 



J ds J ln[(t - s) 2 + {z - r) 2 }dT. 



-I -H 

For the z - component of the gravitational field one has 

l -H+x(s) 



dV(t,z) 
dz 



-H 



An J (t- s) 2 + (z + H - x(s)) 2 

In particular, on the surface z = we obtain the following nonlinear operator 
equation 

i 

ip(x) = £-J K{t, S , x(s))ds - y(t) = 0, (4.1) 
-i 

where 

K(t,s,x(s)) =ln- ( t - s ) 2+H2 



(t-s) 2 + {H-x{s)) 2 ' 

The gravity strength anomaly y(t) = — dV g^ is given and the interface between two 
media (with and without the sources of a gravimetry field) x(s) is to be determined. 

Let ip act between the pair of Hilbert spaces Hi and H 2 - Assume that Hi = 
H 1 ^ 1,1] or Li\— 1,1] and Hi = Li\— 1,1]. The Frechet derivative of this operator is 
the following one 



For any fixed x G {x G Li2[— I, I], x < H — e, e > 0} the kernel 

2(ff-x(a)) 



(t - s) 2 + (# - x(s)) 2 



is a square integrable function on [—1, 1] x [—1, 1], therefore (p'(x) in (4.2) is a compact 
linear operator in L/2[—l,l]. This means that the operators <p'(x) and tp'* (x)ip' (x) 
are not boundedly invertible. So, one can not use classical iterative schemes such 
as Newton or Gauss - Newton in the case of equation (|4.l| ). 

We solve the Cauchy problem ( |2.l| ) for tp given by ( f4.l| ) with some regularization 
function a(t). The problem is numerically solved by means of two finite difference 
methods, namely, Euler's method 

x k+ i = x k +tF (t k ,x k ), x = x(0) (4.3) 

and the Runge - Kutta method 



x k+ i = x k + -F(t k , x k ), (4.4) 



x k+1 = X k + TF(t k+ l,X k+ l), x = x(0). 

Here 

F(t k ,x k ) = -[ip'*(x k )ip'(x k ) + a(t k )I]^ 1 (p'*(x k )ip(x k ) + a(t k )(x k - x )) 
and an equal grid size r > defines the node points, 

t k = fcr, k = 0, 1, .. 

For the successful realization of Continuous Gauss-Newton method an appropriate 
regularization function should be chosen. At the beginning of the process values 
of a(t) should not be very small for the operator ip'* (x(t))ip' (x(t)) + a(t)I to be 
stably invertible and at the same time a(t) should tend to zero sufficiently fast to 



ensure convergence of the function x(t) to the solution of problem (hi). In the 
numerical experiments the functions oto/(fi + t) m , a^e~^ 1 and ao2~ /3t were used. 
The experiments have shown that for all the considered functions the numerical 
solution is evaluated with appropriate accuracy for sufficiently large range of the 
parameters m, ao and (3. The following tables illustrate the dependence of the 
accuracy of the numerical results on parameters for the following data 2 = 1, H = 



2, p = 1, xq — 1. For the numerical tests the function y(t) in (hi) was chosen 
as the solution of the direct problem for the model function x mo d(t) = (1 — t 2 ) 2 . 



The integral in (4.1) was calculated by Simpson's formula with the number of node 
points equals 201 and with a step size equal to 0.01. 

In the tables below and A# are the absolute errors, <je and an are the 



discrepancies ||^>(x(f))||, see (4.1), of the Euler and the Runge-Kutta methods re- 
spectively. The first table shows the dependence of the absolute errors and the 
discrepancies <p(x) on the regularization function. 



Table 1. 











= 0.1, T = 0.1 












a(t) 


N 


A E 




A R 




a (l + t)~ 


2 


38 


0.26 




3.75 • 10-^ 


0.27 




4.09 


• 10 


-2 


a (l + t)- 


4 


25 


0.25 




0.25 


0.11 







11 




a (l + t)- 


6 


43 


0.12 




0.12 


1.75 • 10 


-2 


2.03 


• 10 


-2 


a (l + t)- 


8 


61 


5.30-10 


-2 


2.87- 10" 3 


5.33-10 


-2 


3.58 


• 10 


-3 


a (l+t)- 


10 


79 


2.25-10 


-2 


4.90- 10~ 4 


2.23-10 


-2 


6.22 


• 10 


-4 






127 


1.07-10 


-2 


1.14- 10~ 5 


1.08-10 


-2 


3.03 


• 10 


-5 


aoe" 354 




85 


1.08-10 


-2 


2.84- 10" 4 


1.08-10 


-2 


3.60 


• 10 


-4 



Then it is assumed that the type of a function a(t) is chosen and the dependence 
of the absolute errors and the discrepancies on parameters (3 (Tables 2 and 3) and 
ao (Table 4) is analyzed. 



Table 2. 







a(t) = 


a e 


~P , 


a 


= 0.1, T = 0.1 








/3 


N 


A E 


OE 


Ar 


<?R 


1 


29 


0.26 




7.60 


• 10 


-2 


0.26 


9.10 


10 


-2 


2 


156 


1.09-10 


-2 


4.58 


• 10 


-6 


1.10- 10~ 2 


1.52 


10 


-5 


3 


100 


1.06-10 


-2 


5.78 


• 10 


-5 


1.07- 10~ 2 


9.30 


10 


-5 


4 


73 


1.15 • 10 


-2 


7.54 


• 10 


-4 


1.14- 10~ 2 


1.10 


10 


-3 


5 


56 


1.52 • 10 


-2 


4.50 


• 10 


-3 


1.50-10- 2 


6.12 


10 


-3 


6 


45 


2.08-10 


-2 


1.40 


• 10 


-2 


2.17-10-2 


1.84 


10 


-2 


7 


37 


3.90-10 


-2 


3.06 


• 10 


-2 


3.89-10-2 


3.80 


10 


-2 


8 


31 


0.11 




5.82 


• 10 


-2 


8.59-10-2 


6.52 


10 


-2 


9 


26 


0.18 







10 




0.15 





11 




10 


23 


0.24 







14 




0.20 





16 





Table 3. 







a(t) = 


a e tit , a 


= 0.1, T = 


0.6 








& 


N 


A £ 








1 


4 


0.23 




5.35 • 10" 


-2 


0.26 







12 




2 


20 


4.13- 10 


-2 


1.00-10" 


-3 


1.09-10 


-2 


6.11 


10" 


-5 


3 


16 


8.90-10 


-2 


5.04-10" 


-2 


1.05-10 


-2 


1.73 


10" 


-4 


4 


12 


1.46- 10 


-2 


1.29-10" 


-4 


1.25-10 


-2 


1.79 


10" 


-3 


5 


9 


1.76-10 


-2 


7.14-10" 


-4 


1.84-10 


-2 


1.08 


10" 


-2 


6 


7 


0.16 




3.85 • 10" 


-3 


2.98-10 


-2 


3.22 


10" 


-2 


7 


6 


0.38 




2.58-10" 


-2 


0.12 




5.61 


10" 


-2 


8 


5 


0.83 




1.48 




0.18 




8.60 


10" 


-2 


9 


4 


0.83 




01.48 




0.21 







15 




10 


3 


0.83 




01.48 




0.23 







19 





Table 4. 







a(t) = a a e 


/3 = 3. 


5, r = 0.1 






N 


A E 




Ak 


OR 


10" 3 


71 


1.62 • 10~' 2 


9.77- IO" 4 


1.61 • IO" 2 


1.33- 10" z 


io- 2 


79 


1.33 • IO" 2 


4.55- IO" 4 


2.29- IO" 2 


2.90- IO" 4 


lO^ 1 


85 


1.07- IO" 2 


2.84- IO" 4 


1.08- IO" 2 


3.60- IO" 4 



Analyzing the results of the numerical experiments (a part of them is included 
in the Tables) one concludes the following: 

(i) The Rungc-Kutta method is more stable with respect to changes of the regu- 

larization function a(t) and the step size r, than the Euler method; 

(ii) for all the considered functions suitable parameters can be chosen, however in 

the case when a(t) = u a e~P t the accuracy with which the solution x(t) is 
calculated is higher; 

(iii) an appropriate range of values of the parameter ao is from 0.001 to 0.1, for 
larger values the accuracy is lower, and for smaller values the processes do not 
converge; 

(iv) the range of appropriate values of (3 is large enough: from 2 to 6 for t = 0.6 
and from 2 to 7 for r = 0.1. 



Remark 4.1. The reason why result (i), formulated above, is emphasized can be 
understood if one remembers that problem (4.1) is ill-posed. If a problem is ill- 
posed then the usage of a higher-order accuracy difference scheme (or quadrature 
formula) may lead to less accurate results, as was observed in the literature (see 
e.g. (a), p. 155). On the other hand, if a problem is well-posed, then the usage of 
a higher-order accuracy scheme should lead to more accurate results. 
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